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Abstract 


A 3D simulation tool for solid oxide fuel cells (SOFCs) was described to simulate the mass, momentum and energy conversions in the mono- 
block layers built (MOLB)-type SOFC system. Considering the co-flow and counter-flow cell designs, the temperature distributions, variations of 
reaction species and current densities of the single-unit cell were calculated under the different working conditions. The simulation results show 
that the co-flow case has more uniform temperature and current density distributions. Similar to the planar SOFC, in co-flow case, increasing fuel 
delivery rate or hydrogen mass fraction in the fuel, average temperatures of PEN (positive/electrolyte/negative) and current densities rise, but the 
average temperatures of PEN decrease with increasing the delivery rate of air. In particular, MOLB-type SOFC has some advantages such as: 
higher hydrogen utilizations, lower temperature difference and higher current density. However the current density distributions are less uniform 


in MOLB-type SOFC, which is a disadvantage in this type SOFC. 
© 2007 Published by Elsevier B.V. 
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1. Introduction 


The solid oxide fuel cell (SOFC) is expected to be a promis- 
ing alternative power source for distributed or residential power 
plants because of its higher energy conversion efficiency and 
power density, lower environmental hazards and production 
cost [1-3]. However, the further development of SOFC stacks 
faces the challenges related to maximize the power density and 
to minimize the non-uniform temperature distribution, which 
contributes to the thermal stress in the SOFC components 
[4,5]. 

At this stage of SOFC stack development, design of geometry 
is almost as crucial as material development. In the past, some 
focus has been placed on the designs of new material [6—9] or 
novel geometry [10-13]. Besides planar and tubular SOFCs, 
recently, MOLB-type SOFC also has been drawn keen attention 
and researched to improve the comprehensive performance of 
SOFC. In the MOLB-type SOFC, the corrugated-shaped PEN 
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not only ensures the effective electrochemical reaction area 
is higher than the projected area, but also provides the film 
with the combined function of fuel and airflow paths, mak- 
ing the cell compact and saving the laborious work of channel 
machining. 

In addition to the geometrical designs of SOFC, the work- 
ing conditions, such as delivery rate of gas to the cell system 
(including the fuel and air gas), and hydrogen mass fraction in 
fuel gas also influence the performances of SOFC in a com- 
plicated way. Therefore, it is difficult and uneconomical to 
investigate the influences of these parameters independently by 
testing methods. The computer simulation technique has been 
used to analyze effectively the process of energy conversion in 
the SOFC system. Some modeling of the SOFC during steady 
operation has been constructed to calculate temperature and cur- 
rent density distributions [14—20]. Investigations of planar and 
MOLB-type SOFC operation and performance have predicted 
cell temperature and current density distributions for various 
flow patterns [21-24]. However, only a little work has been 
performed on the influences of the operating conditions on the 
performances of planar SOFC [25,26], or even saying nothing 
of the MOLB-type SOFC. 
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Nomenclature 


Cp specific heat (Jkg~! K7!) 


Ck concentration of component k 

Dxett effect diffusivity of component k (m2 s7!) 

Et internal energy of mixture gas (J) 

Es internal energy of solid (J) 

F Faraday constant (C mol`!) 

i local current density (A cm~?) 

Ik mass source of mixture component k 

j transfer current density (A cm?) 

keff effect thermal conductivity coefficient (W m7!) 
ke thermal conductivity coefficient of gas (W m7!) 
ks thermal conductivity coefficient of solid (W m~ 1) 
K permeability coefficient (m7) 

P pressure of mixture (Pa) 

R gas constant (J mol`! K7!) 

Sk reaction coefficient of component k 

SE energy source (J) 

SM momentum source 

AS entropy change (J mol~!) 

T temperature (K) 

U mixture velocity (m s7 1) 


Greek letters 


E porosity 

ô anode thickness (m) 

Nact,a activation potential at anode (V) 

Nact,c activation potential at cathode (V) 

Left effect viscosity 

p density (kg m7?) 

Oeff effect electrical conductivity (Q27! m7!) 


The objective of present work is to simulate the thermal 
and electrochemical performances of MOLB-type SOFC. A 
computational fluid dynamics (CFD) model tool is demonstrated 
to predict temperature, species mass fraction and current den- 
sity distributions in SOFC system. The results simulated in this 
paper can not only guide the designer in understanding how 
geometrical design affects thermodynamics performances and 
electrochemical characteristics of the SOFC, but also provide 
a more reasonable basis for designing geometry and choosing 
working conditions of the SOFC stacks. 


2. Model development 
2.1. Geometry of MOLB-type SOFC 


As shown in Fig. 1, a typical MOLB-type SOFC [13] is sym- 
metric about the mid-planes of the air and fuel channels. For 
the sake of simplicity in the simulation, the calculation domain 
includes the region between these two mid-planes [21]. The 
thicknesses of anode, cathode, electrolyte and inter-connector 
were 0.5, 0.25, 0.05 and 1.0mm, respectively. 


(a) 


(b) 


Inter-connector 


0.05mm Electrolyte 
0.5mm Anode 


| 0.25mm Cathode 


0.5mm 


Fig. 1. Illustrations of the one cell-stack and the single unit model for MOLB- 
type SOFC. 


2.2. Thermo-fluid model 


In the simulation, the solid and fluid domains were divided 
into some discrete meshes, and in each computational mesh, the 
conservation equations of species, mass, momentum and energy 
were solved using the finite volume method. 


The species conservation equation: 
V(pCxU) = V(Dx,ct#VCx) + Ik, k = H2, O2, H20 (1) 


where J; is the rate of production or the consumption of specie 
k, and given by [27]: 
(Ki) 

QF) 

The diffusion coefficients of hydrogen and oxygen are obtained 
as follow [27]: 


(2) 


I= 4 


T \15 T \15 
Do, =0.181(5) , Dp, =0.753( =) (3) 
The mass conservation equation: 
V -(epU) =0 (4) 


The gas was considered as ideal gas mixtures with the densities 
given by: 


O) ° 
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where m, is the mass fraction of specie k with molecular weight 
Mx. 
The momentum conservation equation: 


. ou p J?” p ðw 
u w 
P ax dy əz 


3u? æ? Ow? 
axe + ay? + a + Sm (6) 


= —EV P + &ett ( 


where Sm is momentum source, and SĮ =0 in the flow chan- 
nels. However, in the porous electrodes, Darcy law with constant 
porosity and permeability is applied to the model and the 
momentum source is obtained by [27]: 


Su = — (=) £U (7) 


where eff is the effective viscosity of the mixture gas and is 
given by [27]: 


X kik 
=>. | 8 
ie 2 Gren) - 


jes (4 + w/w wii) 
[8 + (Mx /M,))| 7 


(9) 


where X; is the molar percent of the specie k, uj and px are 
kinematical viscosities of specie j and k, respectively. 
The energy conservation equation: 


V-(U(prEs + p)) + V (tU) + V: (kerVT) + SE=O (10) 


Heat transfer between the fluid and the solid materials was lim- 
ited to conduction and convection, and the effect of radiation 
was neglected in this calculation because it is very small rela- 
tive to the other kinds of heat transfer. Additionally, the effective 
thermal conductivities of porous electrodes are calculated by the 
Eq. (11) [28]: 


keff = eke + (1 — £)ks (11) 


where kf and ks are thermal conductivities of fluid and solid, 
respectively. Sp is the energy source expressed by the Eq. (12) 
and mainly consists of reaction and Ohmic heats [27]. 


i i i TAS p 
e= + (i) (F) ma) (12) 


Specific heats of gas components are described as the functions 
of temperature [27]: 


Table 1 

Coefficients of the specific heats of gas 

Gas a b c 
Hydrogen 25.8911 —0.8373 2.0138 
Oxygen 29.0856 12.9874 —3.8644 
Water gas 30.3794 9.6212 1.1848 


where a, b and c are relevant coefficients, as listed in Table 1. 
Solid material properties used in this simulation are listed in 
Table 2. 


2.3. Electrochemical model 


The oxidant reduction reaction occurring at the cathode is 
expressed as follows: 


0.50) + 2e7 > O77 (14) 


The oxygen ions transfer through the electrolyte and then into 
the active reaction areas of anode. The electrochemical reaction 
of fuel at the anode is: 


H + 07> > H2O + 2e7 (15) 
So the overall reaction is: 
H2 + 0.502 > H20 (16) 


According to the Faraday law, the reaction rates depend on 
the current density 7 [15]: 


d 
a = 4k (17) 


where df/dt, dO2/dt are the molar consumption rates of fuel and 
oxygen at the anode and the cathode, respectively. 

During the process of energy conversion, when the charge 
transfer reactions at the interfaces of electrolyte and electrode 
are too slow to provide ions at the rate required by the demands 
of current, the activation polarization occurs and is defined by 
the Butler—Volmer equation [28]: 


ae 2F 2F 
t= l0 {exp (-» (=) na = fewa — B) (=) na \ 
(18) 


Eq. (18) is simplified and described by the Tafel empirical for- 
mula [28,29]: 


RT 
= i 19 
Cp =atbx 10T +c x 10°T? (a3) "eeta (sin) na 
Table 2 
Properties of the solid material 
Cell component Density (kg m7?) Effect thermal conductivity Specific heat Porosity (%) Permeability 
(Wm! K!) (Ikg-! K7!) coefficient (m?) 
Inter-connector 7700 13 0.8 - - 
Anode 6200 6.23 0.65 35 1.0E—12 
Cathode 6000 9.6 0.9 35 1.0E—12 
Electrolyte 5560 2.7 0.3 - -= 
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(20) 


RT . RT . 
Nact,c = = (sr) In ioc + (r) lni 
where £ is the transmission coefficient and 6 = 0.5 in this simula- 
tion, Nact,a and Nact,c are the activation potentials at the anode and 
the cathode, respectively. io, and ig, are the exchange current 
densities at the anode and the cathode, respectively. 

A simple semi-empirical formula is used to obtained the 
Nernst voltage: 


RT : 
E = Eo + In ( 720. 
CF) | PH Po; 
0.5 
RT H 
= ink GR in | a (21) 
(2F) PHO 


where Ep is the standard voltage of the cell, PH,0,Po, and Pp, 
are the partial pressures of water gas, oxygen and hydrogen, 
respectively. 


3. Numerical implementation 


The cell unit analyzed in the paper represented a repeating 
unit in the middle of a large stack, and external walls of the 
cell unit were assumed to be adiabatic. Constant temperature, 
delivery rate, and gaseous composition were imposed at the inlet 
boundaries for the fuel and air. 

In the calculations, the modeling tool coupled a thermal-fluid 
model with an electrochemical model. The thermal-fluid model 
was implemented via the commercial CFD simulation code. 
First, the finite-volume Navier-Stokes and transport equations 
were solved to obtain the gas species concentrations and tem- 
peratures at each position in the cell. Then, the information was 
passed to the electrochemical model, which was called via the 
subroutine. Using this solution, the Nernst voltage and the cur- 
rent density distribution were calculated and applied to obtain 
heat source and species source. Finally, gas species concentra- 
tions and temperature distributions were then calculated again 
and provide for the next iteration. The models were coupled time 
after time until convergence of solution was achieved. 


Table 3 
The cell operating conditions and parameters used for simulation 


(a) Temperature [K] 
— 1.063e + 003 


— 1.054e + 003 
— 1.045e + 003 
1.037e + 003 
1.028e + 003 
1.019e + 003 
1.010e + 003 


1.001e + 003 


9.925e + 002 Fuel flow 


9.837e + 002 Air flow 


(b) Temperature [K] Air flow 


1.094e + 003 


1.083e + 003 


1.073e + 003 
1.062e + 003 
1.051e + 003 
1.041e + 003 
1.030e + 003 
1.019e + 003 
1.009e + 003 


9.980e + 002 Fuel flow 


Fig. 2. PEN temperature distributions in the co-flow case (a) and counter-flow 
case and (b) under first working conditions. 


4. Simulation results and discussions 
4.1. Thermo-fluid analysis 


Fig. 2 shows the PEN temperature distributions in the MOLB- 
type SOFC under the first working conditions illustrated in 
Table 3. As comparing Fig. 2(a) with Fig. 2(b), two characteris- 
tic features can be seen. One is that the PEN average temperature 
is 1014K with maximum and minimum temperatures of 1063 


Sample number Fuel 


Air Flow pattern 


Delivery rate (v1) Inlet temperature 


Hydrogen mass 


Delivery rate (v2) Inlet temperature (K) 


(ms~!) (K) fraction (%) (ms!) 
1 0.5 973 0.8 3 873 Counter-flow, Co-flow 
0.5 
2 1.0 973 0.8 3 873 Co-flow 
1.5 
0.8 
3 0.5 973 0.9 3 873 Co-flow 
1 
3 
4 0.5 973 0.8 2 873 Co-flow 
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(a) Hydrogen Mass Fraction 
6.571e - 001 


5.896e - 001 
5.221e - 001 
4.546e - 001 
3.871e - 001 


3.195e - 001 


2.520e - 001 


— 1.845e - 001 


| 


Fuel flow 
Air flow 


= 1.170e - 001 


= 4.952e - 001 


(b) Hydrogen Mass Fraction 


6.480e - 001 A 
Air flow 


5.865e - 001 
5.249e - 001 
4.634e - 001 
4.019e - 001 
3.403e - 001 
2.788e - 001 
— 2.173e - 001 
= 1.557e - 001 


= 9.420e - 002 Fuel flow 


Fig. 3. Hydrogen mass fraction distributions on the interface between 
anode/electrolyte in the co-flow case (a) and counter-flow case (b) under first 
working conditions. 


and 987 K in MOLB-type SOFC in co-flow case (see Fig. 2(a)). 
But, in counter-flow case, average temperature is 996K with 
maximum and minimum temperatures of 1088 and 923 K (see 
Fig. 2(b)). So the co-flow case has the more uniform tempera- 
ture distribution and smaller temperature difference (AT) from 
air inlet to outlet. A further difference is that the temperature 
of PEN increases uniformly along the direction of fuel flow, 
and is highest near the fuel outlet in co-flow case. However, the 
temperature of PEN rises rapidly, reaching a maximum near the 
fuel inlet, and then gradually drops for counter-flow case. This 
is due to the offsetting effects of air near the inlet, at its coolest, 


being aligned with the fuel inlet. As a consequence, the temper- 
ature gradient is smaller in co-flow case, which must result in 
the smaller thermal stress in the PEN although the thermal stress 
distributions were not researched in the paper. 

Fig. 3(a and b) illustrates the hydrogen mass fraction distribu- 
tions in the co-flow and counter-flow cases, respectively. Along 
the direction of fuel flow, the mass fractions of hydrogen on the 
interface between the electrolyte and anode decrease due to the 
electrochemical reaction. In particular, of the two flow case, the 
hydrogen mass fraction near the fuel outlet is less in co-flow 
case, so the more fuel is consumed and fuel utilization is higher 
accordingly. 

On the basis of above analysis about the temperature distri- 
bution, it was found that the co-flow case is advantageous to 
improve the fuel utilization and mitigate the steep temperature 
gradient, and hence to reduce the internal stresses. For the co- 
flow case, other working conditions influenced the performances 
of SOFC were also studied. Firstly, we characterized the effects 
of the fuel gas on current density and temperature distributions. 
Fig. 4 shows the current density distributions of PEN. Changing 
the delivery rate of fuel with 0.5, 1.0 and 1.5 m s7}, respectively, 
the average current density gradually increases (see Fig. 4(a—c)), 
but the current density is less uniform. This indicates that it is 
effective for the improvement of the electrical performance to 
increase the delivery rate of fuel. Fig. 5(a) shows the temperature 
distributions of the mid-plane in the X-direction. With the incre- 
ment of the delivery rate of fuel, the average temperature and 
temperature difference (AT) from fuel inlet to outlet of PEN rise, 
so the temperature gradient also rises, which may cause thermal 
stress. 

Then, we focus on the influence of hydrogen proportion in 
fuel on the temperature and current density distributions. Chang- 
ing the hydrogen mass fraction in the fuel gas with 0.8, 0.9 
and 1.0, reaction rate rise and more hydrogen is consumed, and 
then the more reaction heat is accumulated. As a result, the 
temperature gradient and current density also rise, as shown in 
Fig. 5(b) and Fig. 6(a—c). Furthermore, as comparing Fig. 4 with 
Fig. 6, with higher delivery rate of hydrogen, current density 
distributions were less uniform, but with larger mass fraction 
of hydrogen in fuel gas, current density distributions are more 
uniform. 

In order to decrease the temperature gradient, it is effective 
to increase delivery rate of air. In the SOFC system, air not 
only provides oxygen ions but also has the cooling function. 


(a) Current Density [A cm*-2] (b) Current Density [A cm*-2] (c) Current Density [A cm4-2] 
-1.192e + 000 -1.241e + 000 -1.272e + 000 : 
-1.126e + 000 -1.174e + 000 -1.208e + 000 
-1.060e + 000 -1.107e + 000 -1.144e + 000 
-9.943e - 001 1.040e + 000 -1.080e + 000 

9.283e - 001 9.730e - 001 1.016e + 000 

8.623e - 001 9.061e - 001 9.517e - 001 

7.963e - 001 8.391e - 001 8.877e - 001 

7.303e - 001 7.721e - 001 8.236e - 001 

6.644e - 001 Fuel flow $ 7-051e - 001 Fuel flow 7.595e - 001 Fuel flow 
5.984e - 001 Air flow 6.381e - 001 Air flow 6.955e - 001 Air flow 


Fig. 4. Current density distributions on the interface of anode/electrolyte for the MOLB-type SOFC in the co-flow case under second working conditions. 
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(b) 1180] —*—H2.mf1.0 cyte 
1160 —e—H,.mf=0.9 ~ 
—s—H,.mf=0.8 


0.2 0.4 0.6 0.8 
Distance from fuel inlet (dm) 
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Fig. 5. PEN temperature distributions for MOLB-type SOFC under the second working conditions (a) and third working conditions (b). 


(a) Current Density [A cm^-2] 


-1.192e + 000 - 1.115e + 000 

-1.126e + 000 ~ 1.073e + 000 

~ 1.060e + 000 m 1.031e + 000 
9.943e - 001 9.882e - 001 
9.283e - 001 9.458e - 001 
8.623e - 001 9.034e - 001 
7.963e - 001 8.610e - 001 
7.303e - 001 8.186e - 001 
6.644e - 001 Fuel flow 7.762e - 001 
5.984e - 001 Air flow 7.338e - 001 


(b) Current Density [A cm^-2] 


(c) Current Density [A cm^-2] 
-1.210e + 000 


-1.170e + 000 
-1.131e + 000 
1.092e + 000 
1.053e + 000 
1.013e - 000 
9.739e - 001 
9.346e - 001 
8.953e - 001 
8.560e - 001 


j 


Fuel flow 
Air flow 


f 


Fuel flow 
Air flow 


Fig. 6. Current density distributions on the interface of anode/electrolyte for MOLB-type SOFC under third working conditions. 


Increasing the delivery rate of air, the steep temperature gradient 
is mitigated in MOLB-type SOFCs, as shown in Fig. 7. This is 
because more reaction heat is absorbed and then released by the 
air with the higher delivery rate, although the air utilization is 
dropped. 


4.2. Advantage of MOLB-type SOFC compared with planar 
SOFC 


On the basis of the preceding works about the simulating the 
planar SOFC system [26], the general trends of the temperature 
distributions in the planar and MOLB-type SOFCs are similar 
in essential. Although the average temperatures in MOLB-type 
SOFC are higher than those in the planar, the temperature dif- 
ference (AT) is lower in MOLB-type SOFC. This is because the 
electrochemical reaction area of MOLB-type SOFC is larger and 


(a) Temperature [K] 
- 1.084e + 003 


~ 1.074e + 003 
~ 1.063e + 003 
1.053e + 003 
1.042e + 003 
1.032e + 003 
1.022e + 003 
1.011e + 003 
1.001e + 003 
9.904e + 002 


} 
Fuel flow 
Air flow 


the more fuel is consumed in the active areas, as result, the more 
hydrogen is consumed and more reaction heat is accumulated. 
However, due to the effect of corner formed by the inclined plane 
and the upper plane, the temperature and hydrogen mass fraction 
distributions in MOLB-type SOFC are less uniform than those 
in the planar. 

As comparing the current density distribution in planar SOFC 
with that in the MOLB-type SOFC, some main characteristic 
features can be seen. The first feature is that, in the two kinds 
of SOFCs, the current densities gradually decrease along the 
fuel flow direction because of the dropping of hydrogen mass 
fraction. And the second is that, average current densities are 
higher than those in the planar SOFC due to the larger reaction 
areas in MOLB-type one, but the current density distributions 
are less uniform in the MOLB-type SOFC. This is because, the 
current densities not only decrease along the fuel flow direction 


(b) Temperature [K] 
- 1.104e + 003 


~ 1.093e + 003 
= 1.083e + 003 
1.072e + 003 
1.061e + 003 
1.051e + 003 
1.040e + 003 
1.029e + 003 
1.019e + 003 
1.008e + 003 


Fuel flow 
Air flow 


Fig. 7. PEN temperature distributions for MOLB-type SOFC under forth working conditions. 
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but also change along the incline plane and reach the maximum 
on the top plane due to the gathering current density function 
of the inter-connector. The last characteristic is that, with the 
increase in the delivery rate and mass fraction of hydrogen, the 
increment of the temperatures and current densities are less than 
those in planar SOFC. 


5. Conclusions 


The exchange current density was used to couple the thermo- 
fluid model with electrochemical model. The performances of 
MOLB-type SOFCs in stead state were calculated and com- 
pared with the performances of the planar SOFC. Some main 
conclusions were made as follows: 

The PEN temperature distribution in the MOLB-type one is 
more uniform in co-flow case than that in counter-flow case, in 
particular, the temperature gradient is lower and current density 
is higher than those in the planar SOFC. 

For co-flow case, with increasing the delivery rate of fuel 
gas or hydrogen mass fraction in the fuel, temperature gradients 
and average current densities rise in the MOLB-type SOFC. But 
increasing delivery rate of air, steep temperature gradients are 
mitigated. In particular, with the same increasing in delivery rate 
of fuel and hydrogen mass fraction, the increment of temperature 
distributions in MOLB-type SOFC are less. The current density 
distributions are less uniform with higher delivery rate of fuel, 
and are more uniform with larger mass fraction of hydrogen in 
fuel. 
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